Primeiro, programo os métodos numéricos. Crio um método Runge-Kutta genérico, que aceita uma matriz de butcher como parâmetro, e depois ponho rk4 como uma aplicação parcial deste método com uma matriz especÃfica.
library(tidyverse)
library(gganimate)
library(udunits2)
rungeKutta = function(M, f, y0, ts) {
y = matrix(0, nrow = length(ts), ncol = length(y0))
colnames(y) = if (is.null(names(y0))) paste0("y", 1:length(y0)) else names(y0)
y[1,] = y0
# butcher tableau
c = M[1:(nrow(M)-1), 1]
A = M[1:(nrow(M)-1), 2:ncol(M)]
b = M[nrow(M), 2:ncol(M)]
for (i in 1:(length(ts) - 1)) {
h = ts[i + 1] - ts[i]
K = matrix(0, nrow = nrow(A), ncol = length(y0))
K[1,] = f(ts[i], y[i,])
for (j in 2:nrow(K)) {
K[j,] = f(ts[i] + h*c[j], y[i,] + h * as.vector(A[j, 1:(j-1)] %*% K[1:j-1,]))
}
y[i + 1,] = y[i,] + h*b %*% K
}
y
}
rk4 = function(f, y0, ts) {
M = matrix(c(0, 0, 0, 0, 0,
.5, .5, 0, 0, 0,
.5, 0, .5, 0, 0,
1, 0, 0, 1, 0,
0, 1/6, 1/3, 1/3, 1/6), byrow = TRUE, nrow = 5)
rungeKutta(M, f, y0, ts)
}
Precisamos de dados. Coletei os parâmetros padrão de gravitação da wikipedia e os ephemeris da Horizons - NASA. Precisamos fazer conversões de unidades, pois os parâmetros de gravitação estão em \(m^3s^{-2}\) e os ephemeris, em unidades astronômicas. Uso como unidade de distância \(10^6 km\) por conveniência.
init = read_csv("../data/ephemerides-2018-09-26.csv") %>%
mutate(sgp = ud.convert(sgp, "m^3/s^2", "(1e6*km)^3/(days^2)")) %>%
mutate_at(vars(x, y, z), ~ ud.convert(., "au", "1e6km")) %>%
mutate_at(vars(vx, vy, vz), ~ ud.convert(., "au/days", "1e6km/days"))
init
## # A tibble: 10 x 8
## Body sgp x y z vx vy vz
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Sun 9.93e+2 -1.71e-3 1.08e0 -0.0114 -0.00113 4.13e-4 2.81e-5
## 2 Mercu… 1.64e-4 -5.89e+1 -1.63e1 3.98 0.317 -3.86e+0 3.44e-1
## 3 Venus 2.43e-3 1.05e+2 -2.86e1 -6.45 0.811 2.90e+0 -7.06e-3
## 4 Earth 2.98e-3 1.50e+2 7.90e0 -0.0119 -0.160 2.56e+0 -5.21e-5
## 5 Mars 3.20e-4 1.97e+2 -6.21e1 -6.16 0.718 2.17e+0 2.78e-2
## 6 Jupit… 9.48e-1 -4.16e+2 -6.88e2 12.1 0.953 -5.30e-1 -1.91e-2
## 7 Saturn 2.83e-1 2.18e+2 -1.49e3 17.2 0.779 1.18e-1 -3.32e-2
## 8 Uranus 4.32e-2 2.57e+3 1.49e3 -27.8 -0.299 4.82e-1 5.65e-3
## 9 Neptu… 5.10e-2 4.32e+3 -1.16e3 -75.7 0.119 4.56e-1 -1.21e-2
## 10 Pluto 6.51e-6 1.74e+3 -4.73e3 4.05 0.450 6.12e-2 -1.38e-1
A função que passamos para o método RK4 tem como parâmetro o tempo e o estado, que é um vetor com 60 componentes: cada corpo tem parâmetros x, y, z e vx, vy, vz. Ou seja, estamos resolvendo um sistema com 60 equações diferenciais.
Não há muita novidade. A equação é:
\[\frac{d²x_i}{dt²} = \sum_{j \neq i} \frac{\mu_j(x_j-x_i)}{\mid x_j-x_i \mid^3}\]
f = function(t, state) {
n = length(state)/6
x = state[1:n]
y = state[(n+1):(2*n)]
z = state[(2*n+1):(3*n)]
vx = state[(3*n+1):(4*n)]
vy = state[(4*n+1):(5*n)]
vz = state[(5*n+1):(6*n)]
ax = ay = az = double(n)
for (i in 1:n) {
dx = x[-i] - x[i]
dy = y[-i] - y[i]
dz = z[-i] - z[i]
r2 = dx^2 + dy^2 + dz^2
ax[i] = init$sgp[-i] * dx / r2^1.5
ay[i] = init$sgp[-i] * dy / r2^1.5
az[i] = init$sgp[-i] * dz / r2^1.5
}
c(vx, vy, vz, ax, ay, az)
}
Por conveniência, crio um wrapper que gera os resultados, para uma sequência \(t_n\), já no formato adequado para os plots.
make_tb = function(ts) {
state0 = c(init$x, init$y, init$z, init$vx, init$vy, init$vz)
res = rk4(f, state0, ts)
tb = map_dfr(1:length(ts), function(i) {
state = init
n = nrow(state)
state$x = res[i, 1:n]
state$y = res[i, (n+1):(2*n)]
state$z = res[i, (2*n+1):(3*n)]
state$vx = res[i, (3*n+1):(4*n)]
state$vy = res[i, (4*n+1):(5*n)]
state$vz = res[i, (5*n+1):(6*n)]
state$t = ts[i]
state
})
}
Gráficos para os planetas mais próximos ao sol (até marte):
tb = make_tb(seq(0, 365*3, by = 2))
write_rds(tb, "published/short_run.rds")
tb %>%
ggplot(aes(x = x, y = y, color = Body)) + geom_point() + theme_dark() +
xlim(-250, 250) + ylim(-250, 250) + xlab("") + ylab("") +
theme(
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
) +
transition_time(t) + ease_aes('linear') +
labs(title = "Date = {as.Date('2016-09-26') + as.difftime(frame_time, units = 'days')}") +
shadow_wake(1, size = .25, alpha = .25)
anim_save("solar_in_xy.gif", path = "published")
tb = make_tb(seq(0, 365*3, by = 2))
write_rds(tb, "published/short_run.rds")
tb %>%
ggplot(aes(x = x, y = z, color = Body)) + geom_point() + theme_dark() +
xlim(-250, 250) + ylim(-50, 50) + xlab("") + ylab("") +
theme(
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
) +
transition_time(t) + ease_aes('linear') +
labs(title = "Date = {as.Date('2016-09-26') + as.difftime(frame_time, units = 'days')}") +
shadow_wake(1, size = .25, alpha = .25)
anim_save("solar_in_xz.gif", path = "published")
tb %>%
ggplot(aes(x = y, y = z, color = Body)) + geom_point() + theme_dark() +
xlim(-250, 250) + ylim(-50, 50) + xlab("") + ylab("") +
theme(
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
) +
transition_time(t) + ease_aes('linear') +
labs(title = "Date = {as.Date('2016-09-26') + as.difftime(frame_time, units = 'days')}") +
shadow_wake(1, size = .25, alpha = .25)
anim_save("solar_in_yz.gif", path = "published")
yz
Mesma coisa, mais com zoom out. Como plot apenas uma vez a cada 320 dias, já que os perÃodos agora são muito mais longos.
tb = make_tb(seq(0, 365*260, by = 2))
write_rds(tb, "published/long_run.rds")
tb %>%
filter(t %% 160 == 0) %>%
ggplot(aes(x = x, y = y, color = Body)) + geom_point() + theme_dark() +
xlim(-5000, 7000) + ylim(-5000, 7000) + xlab("") + ylab("") +
theme(
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
) +
transition_time(t) + ease_aes('linear') +
labs(title = "Date = {as.Date('2016-09-26') + as.difftime(frame_time, units = 'days')}") +
shadow_wake(1, size = .25, alpha = .25)
anim_save("solar_out_xy.gif", path = "published")
tb %>%
filter(t %% 160 == 0) %>%
ggplot(aes(x = x, y = z, color = Body)) + geom_point() + theme_dark() +
xlim(-5000, 7000) + ylim(-2200, 1400) + xlab("") + ylab("") +
theme(
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
) +
transition_time(t) + ease_aes('linear') +
labs(title = "Date = {as.Date('2016-09-26') + as.difftime(frame_time, units = 'days')}") +
shadow_wake(1, size = .25, alpha = .25)
anim_save("solar_out_xz.gif", path = "published")
tb %>%
filter(t %% 160 == 0) %>%
ggplot(aes(x = y, y = z, color = Body)) + geom_point() + theme_dark() +
xlim(-5000, 7000) + ylim(-2200, 1400) + xlab("") + ylab("") +
theme(
axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
) +
transition_time(t) + ease_aes('linear') +
labs(title = "Date = {as.Date('2016-09-26') + as.difftime(frame_time, units = 'days')}") +
shadow_wake(1, size = .25, alpha = .25)
anim_save("solar_out_yz.gif", path = "published")
Aqui, apresento os resultados utilizando integração verlet. Aqui, a implementação do método.
verlet = function(f, y0, dy0, ts) {
dt = ts[2] - ts[1]
y = matrix(0, nrow = length(ts), ncol = length(y0))
y[1,] = y0
y[2,] = y0 + dy0*dt + f(y0)*dt*dt/2
for (i in 3:length(ts)) {
y[i,] = 2*y[i-1,] - y[i-2,] + f(y[i-1,])*dt*dt
}
y
}
fv = function(state) {
n = length(state)/3
x = state[1:n]
y = state[(n+1):(2*n)]
z = state[(2*n+1):(3*n)]
ax = ay = az = double(n)
for (i in 1:n) {
dx = x[-i] - x[i]
dy = y[-i] - y[i]
dz = z[-i] - z[i]
r2 = dx^2 + dy^2 + dz^2
ax[i] = init$sgp[-i] * dx / r2^1.5
ay[i] = init$sgp[-i] * dy / r2^1.5
az[i] = init$sgp[-i] * dz / r2^1.5
}
c(ax, ay, az)
}
make_tbv = function(ts) {
state0 = c(init$x, init$y, init$z)
velocity0 = c(init$vx, init$vy, init$vz)
res = verlet(fv, state0, velocity0, ts)
map_dfr(1:length(ts), function(i) {
state = init
n = nrow(state)
state$x = res[i, 1:n]
state$y = res[i, (n+1):(2*n)]
state$z = res[i, (2*n+1):(3*n)]
state$t = ts[i]
state
})
}
O código para gerar os plots é, mutatis mutandis, o mesmo. Vamos para os plots: